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Abstract 

We present an extension of the Karman-Howarth theorem to the 
Lagrangian averaged magnetohydrodynamic (LAMHD—a) equations. 
The scaling laws resulting as a corollary of this theorem are studied in 
numerical simulations, as well as the scaling of the longitudinal struc- 
ture function exponents indicative of intermittency. Numerical simula- 
tions for a magnetic Prandtl number equal to unity are presented both 
for freely decaying and for forced two dimensional MHD turbulence, 
solving directly the MHD equations, and employing the LAMHD—a 
equations at 1/2 and 1/4 resolution. Linear scaling of the third-order 
structure function with length is observed. The LAMHD—a equa- 
tions also capture the anomalous scaling of the longitudinal structure 
function exponents up to order 8. 
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1 Introduction 



Turbulent flows are ubiquitous in Nature and their very complexity make 
them hard to understand, let alone predict. One reason is that turbulence 
comes in intermittent "gusts", such as fronts. These gusts are associated 
with their non-Gaussian statistics, which allows strong events occurring in 
the "fat tails" of the statistical distribution of, say, velocity gradients and 
other physical variables. From the theoretical viewpoint, its property of 
intermittency is the sine qua non of turbulence. Intermittency may also 
trigger large scale effects, such as the apparently random reversals of the 
Earth's magnetic field £Q. It plays a role in other large scale events in the 
solar dynamo, and in the atmosphere of the earth. Although intermittency is 
believed to take place at small scales, strong events can affect the dynamics 
of the large scales, specially in systems close to criticality. For instance, 
Ref. j2] shows that local fluctuations of the kinetic helicity can explain phase 
and amplitude variations of the 22-years solar cycle. Intermittency is also 
one possible explanation for the occurrence of the Maunder-like minima of 
solar activity Finally, intermittency is known to affect the transport of 
momentum in atmospheric surface layers [3]. 

Intermittency is a highly spatially and temporally localized phenomenon, 
which thus requires high-resolution instrumentation, be it in the laboratory, 
in atmospheric and geophysical flows or in numerical simulations. In the 
latter case, the lack of adequate computer power (from the standpoint of the 
geo-physicist) implies that modeling of the unresolved small scales must take 
place. However, intermittency in general is not included explicitly in Large 
Eddy Simulations (LES) models of turbulent processes. 

The intermittency of the subgrid-scale (SGS) energy dissipation as is 
usually defined in LES has been studied filtering 3D direct numerical simu- 
lations (DNS) of non- conductive fluids [S], with a Taylor Reynolds number 
R\ = 150; such a Reynolds number is somewhat insufficient for a determi- 
nation of third-order scaling. For the SGS dissipation, they find that the 
Smagorinsky model, the volume-averaged dynamic model, and the similar- 
ity model perform fairly well (e.g. the error in the exponent for p = 7 is 
less than 7%). On the other hand, the constant eddy- viscosity and spectral 
eddy-viscosity models underestimate intermittency beyond p = 4 (compared 
to DNS) while the local and clipped dynamic Smagorinsky models strongly 
overestimate the intermittency beyond p = 4. In models not capturing in- 
termittency properly, the question arises whether the overall statistics of the 
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flow at large scales would be affected by the absence of intermittency; and if 
so, how intermittency should be incorporated. 

There are many models of turbulent flows (see for example the recent 
review in j^J) and there are many turbulent flows for which adequate testing 
of such models is in order. Magnetohydrodynamics, i.e. the coupling of a 
velocity and magnetic field at sub-luminal velocities so that the displacement 
current can be neglected, presents an interesting property, namely that in two 
space dimensions (2D), intermittency occurs as well as in three dimensions 
(3D). This is in contrast with the 2D neutral fluid case for which the conser- 
vation of vorticity leads to an inverse energy cascade to the large scales; in 
the presence of a magnetic field, this conservation is broken by the Lorentz 
force. Since, from a numerical standpoint, much higher Reynolds numbers 
can be achieved in 2D, an intermittent flow can be reached in 2D-MHD with 
adequate scale separation between the energy-containing range, the inertial 
range and the dissipation range. Our 2D MHD tests are able to exhibit a sub- 
stantially larger Reynolds number (up to R\ ~ 1500) than the values listed 
for previous studies. This fact provides an ideal testing ground for models 
of turbulent MHD flows, such as they occur in geophysics and astrophysics: 
magnetic fields are observed in detail in the Earth and Sun environments, 
and are known to be dynamically important as well for the solar-terrestrial 
interactions (the so-called space weather), in the interstellar medium and in 
galaxies. 

Modeling of MHD flows is still under development (see jjj). Most LES 
for hydrodynamic turbulence are based upon self-similarity or universality, 
in that they assume a known power law of the energy spectrum. For MHD, 
the kinetic energy is not a conserved quantity, and this poses a problem for 
the extension of such techniques to the case of MHD. Additional difficulties 
arise from the fact that MHD has several regimes depending on the relative 
strengths of the magnetic and velocity fields, their degree of alignment, and 
whether mechanical or magnetic energy is injected into the flow. However, 
some LES have been developed for particular cases. There exists LES for 
MHD turbulence with some degree of alignment between the velocity and 
magnetic fields [S], dissipative LES which does not model the interactions 
between the two fields jH], and LES for low magnetic Reynolds number [TU] . 
We have recently tested one model which may be more generally applica- 
ble, the Lagrangian averaged magneto-hydrodynamics alpha (LAMHD—a) 
model, both in 2D [TT] and in 3D ^2] and it has been used to examine the 
onset of the dynamo instability when the magnetic Prandtl number (the ra- 
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tio of viscosity v to magnetic diffusivity 77) is small as occurs in liquid 
metals in the laboratory, in the liquid core of the Earth or in the solar con- 
vection zone. In this context, because of the importance of intermittency as 
a fundamental, or even defining property of turbulence, we seek to determine 
to what extent LAMHD—a exhibits intermittency. 

Intermittency is believed to be associated only with a forward cascade of 
energy; that is, the cascade of energy from larger scales to smaller scales, or, 
equivalently, from low wave numbers to high wave numbers. As previously 
mentioned, in determining the extent to which the LAMHD—a model ex- 
hibits intermittency, we shall take advantage of the forward cascade of energy 
which occurs in two-dimensional MHD. We shall first investigate the spectral 
scaling laws for this situation in Sec. 01 then its Karman-Howarth theorem 
in Sec. EJ Section 0] also discusses the modifications of the Karman-Howarth 
theorem for MHD which arise due to the presence of the length scale alpha (a) 
in the LAMHD—a model. The length scale alpha modifies the nonlinearity 
in the motion equation, and one must estimate its observable physical effects. 
In particular, introduction of the length scale alpha modifies the LAMHD—a 
energy spectrum for ka > 1. Section discusses energy conservation for the 
LAMHD—a model and investigates its inviscid energy dissipation anomaly, 
which arises from its scaling laws and its Karman-Howarth theorem. Finally 
in Sec. El we investigate the effects of introducing the length scale alpha on 
the intermittency of the LAMHD—a model solutions for decaying and forced 
turbulence in two dimensions. Just as for Navier-Stokes turbulence, these ef- 
fects emerge in numerical simulations as a scaling anomaly in the higher 
order structure functions of the LAMHD—a model. Section [7| summarizes 
these results. 

2 Background of the LAMHD—a model 
2.1 The LAMHD—a model equations 

The Lagrangian-averaged magnetohydrodynamic alpha, or LAMHD—a model 
was derived by Lagrangian averaging ordinary MHD along particle trajec- 
tories ^3]. Specifically, the LAMHD—a equations arise from Lagrangian- 
averaging Hamilton's principle for incompressible ideal MHD, after using a 
form of Taylor's hypothesis of frozen-in turbulent fluctuations in the Euler- 
Poincare equation for barotropic MHD from Ref. ^3]. When Navier-Stokes 
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viscosity v and diffusivity r\ are included in the standard fashion, the equa- 
tions for the LAMHD—a model emerge as, 

d t u + u s -Vu-B s -VB + (Vu s ) T ■ u + (VB) T ■ B s + Vtt = uAu, (1) 
d t B s + u s ■ VB S - B s • Vu s = 77AB, (2) 
div u s = 0, div B s = . (3) 

In these equations, subscript s denotes the smoothing obtained by inverting 
the Helmholtz relations, 

u = (1 - a 2 A)u s , B = (1 - a 2 M A)B s . (4) 

with Dirichlet boundary conditions, 

u s = and B s = on the boundary . 

The modified total pressure n in the motion equation for the LAMHD—a 
model is defined by 

1 la 2 

7T+-|B s | 2 =p--| Us | 2 -y|Vu s | 2 , (5) 

where p is the mechanical pressure. In these equations, a and % are two 
constant parameters: a characterizes the correlation length between the in- 
stantaneous Lagrangian fluid trajectory and its mean (time average); while 
qm is its magnetic counterpart. These two parameters need not be equal, ab 
initio. The traditional MHD system is obtained by setting both a = and 
a.M — 0. Likewise, the LANS— a incompressible fluid turbulence model is 
obtained by setting B s = 0. These equations may also be obtained through 
a filtering approach, as proposed in the fluid case in [TH] . 



2.2 Energy, momentum, circulation and linkages 

The LAMHD—a system of equations (001) possesses the standard properties 
of a normal fluid theory. For example, the LAMHD—a system monotonically 
dissipates the positive energy 



£= 2 



l - J (u s -u + B s -B)rf 3 x, (6) 
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according to 

d JL = - v ^ (|Vu s | 2 + a 2 |Au s | 2 )d 3 x- rj J (| VB| 2 + a 2 \ AB| 2 ) d 3 x . (7) 

In addition, the LAMHD—a motion equation (JJ may be expressed in con- 
servative form as 

with stress tensor 

T| = (uiui - a 2 u sk , lU ^ - (B t Bi - a 2 B sk>i B^ 

+ ^'(p-^|B s | 2 + B s .B-^|VB s | 2 ). (9) 

Thus, the two velocities appearing in the LAMHD—a model may be inter- 
preted as fluid parcel velocity u s and momentum per unit mass u. 

The Kelvin circulation theorem for the incompressible LAMHD—a mo- 
tion equation JI} involves both of these velocities, 

u • dx. = (f (J x B s + z/Au) • dx , (10) 



dt Jc(u a ) Jc(u s ) 

where J = curl B. Hence, the J x B s force and viscous force can each generate 
circulation of u around material loops moving with smoothed velocity u s . 
This results by Stokes theorem in vorticity dynamics for uj = curlu in the 
form 

duj 

— + u s ■ Voj - uj ■ Vu s = B s VJ — J VB S + vAu: . (11) 
at 

The linkages of the smooth B-field B s with itself and with the vorticity u) are 
given respectively by the helicity J A s ■ B s d 3 x and cross helicity J u-~B s d 3 x. 
The densities for these linkages satisfy 

-|(A S • B s ) + div ((A, ■ B s ) u s ) = 7] (A ■ B s + A s ■ B) , (12) 



in which B s = curl A s and B = curl A, and, cf. (J3J), 

d_ 
dt 



|^(u-B s )+ div((u-B s )u s + (7r + ^|B s | 2 )B s ) = z/B s Au + r/u-AB. (13) 
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Thus, resistivity affects the helicity, while both resistivity and viscosity affect 
the cross helicity, and these linkages are both preserved by LAMHD—a in 
the ideal case. Of course, these properties of energy, momentum, circulation 
and linkages for the LAMHD—a model all reduce to characteristics of normal 
MHD, when a 2 -> 0. 



2.3 Recasting LAMHD—a as an LES turbulence model 

The LAMHD—a model modifies the motion equation for ordinary MHD. By 
a short sequence of manipulations, we may recast the LAMHD—a motion 
equation into a form which is reminiscent of an LES turbulence model. We 
begin with the following commutation relation, 

[p ■ V, (1 - a 2 A)] q = a 2 div (Vq • Vp + Vq • Vp T ) - a 2 (V(div p) • V) q(14) 
which holds for any vectors p and q. Two other useful vector identities are, 
(Vu s ) T . u - V (||u s | 2 + % | Vu s | 2 ) = -a 2 div (VuJ • Vu s ), (15) 



:vb) t • B s 



V (i (1 - a 2 A) \B S \ 2 + £| VB S | 2 ) = « 2 div (VBj • VB S ) (16) 



where a = au was assumed. Consequently, the motion equation in the 
incompressible LAMHD—a model may be rewritten equivalently in "LES 
form" as, 

(1 - a 2 A) (d t u s + u s ■ Vu s - B s ■ VB S + Vp s - uAu s ) = -a 2 div r, (17) 

where the divergence of the "stress tensor" r is given by 

div r = div (Vu s ■ Vu s + Vu s ■ VuJ - VuJ ■ Vu s ) 

- div (VB S ■ VB S + VB S • VBj - VBj • VB S ) , (18) 

and gradient terms have been absorbed into the modified total pressure, 
denoted by tt, which is given by 

* = P ~ \\»s\ 2 - y |Vu s | 2 - \ (1 - a 2 A) \B S \ 2 - y |VB S | 2 , (19) 
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where p = (1 — a 2 A)p s . By using the following identity for divergenceless 
vectors div u = 

div (Vu T • Vu T ) = V- tr (Vu ■ Vu) , (20) 

we may rewrite the added stress in (JTSjl equivalently, as 

divr = 4div(S-fi-S- J) + Vn (21) 

with new notation 

S = i(Vu s + VuJ), 0=^(Vu s -VuJ) 
£ = ^(VB s + VBj), J=i(VB s -VBj), 

and additional pressure 

II = - tr (Vu s ■ Vu s ) - - tr (VB S • VB S ) . 

The (non-symmetric) stress tensor r given by 

r = 4(S-fi-E- J) + Idn (22) 

which emerges from these manipulations casts the LAMHD—a model into a 
form reminiscent of an LES turbulence model. We shall find these expres- 
sions convenient below in introducing the analog of Elsasser variables for the 
LAMHD—a model. 

3 Inert ial range scaling laws in forward en- 
ergy cascade 

In two space dimensions, the conservation of vorticity in the neutral (B = 
case) leads to an inverse energy cascade; however, the Lorentz force breaks 
this conservation and, in MHD, energy is found to be mostly transferred to 
the small scales both in 2D and in 3D. Several measurements, starting with 
satellite data in the solar wind and continuing more recently with direct nu- 
merical simulations both in two dimensions and three dimensions, indicate 
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that the energy spectrum of a turbulent MHD flow follows a law that is 
barely distinguishable from a neutral fluid, with E(k) ~ k" 1 ' 70 . Differences 
do occur when one examines higher order structure functions: the most inter- 
mittent case (almost comparable in magnitude to that of the passive scalar) 
is the two-dimensional MHD fluid; three-dimensional MHD appears less in- 
termittent than the 2D MHD case ^7] , and the 3D neutral fluid is the lesser 
intermittent of the three. Intermittency has been observed in the Solar Wind 
as well. None of the data at the second order level is in agreement with the 
phenomenologies developed by Kolmogorov jTHj and leading to an energy 
spectrum E{k) ~ k~ 5 ^ 3 (heretofore the K41 model) or by Iroshnikov and by 
Kraichnan ^U] (heretofore, the IK model) and leading to a shallower spec- 
trum, viz. E(K) ~ k~ 3 / 2 in the absence of significant velocity-magnetic field 
correlations. These two types of phenomenology differ by the taking into 
account in the latter case of the non-local interactions (in Fourier space) em- 
anating from the propagation of Alfven waves; it is worth mentioning here 
that the IK model also agrees with the isotropic limit of the weak turbulence 
theory for incompressible MHD [20] • Note also that a model of intermittency 
for MHD flows [21] does recover the intermittency as measured in direct nu- 
merical simulations both in 2D [21] and in 3D J7|, but such models depend 
on two adjustable parameters and thus do not necessarily have a predictive 
power. 

These anomalous scaling laws are not fully understood but, for the neutral 
case, there is an exact law at third order with which the K41 phenomenology 
is compatible. In MHD, the exact law is more complex in its structure since it 
involves third-order cross-correlations between the velocity and the magnetic 
field [22] whereas the phenomenologies evoked above refer to single-variable 
moments. In that instance, it is worth asking what is the equivalent exact 
law in the context of the LAMHD—a model, a task developed in the next 
section. The question also arises as to what is the spectrum of energy beyond 
the alpha cut-off scales (which we take equal here, although different choices 
can be made, see e.g. |12[ 113]). The answer should be guided by what is the 
pseudo-invariant in the small scales, beyond a. In the neutral fluid case, it 
is the enstrophy (u 2 ) and in MHD this becomes (u 2 ) + (j 2 ), where j is the 
current density. A Kolmogorov-like dimensional reasoning (see [23] for the 
LANS— a case) taking into account this pseudo-invariance law will lead to a 
k~ 3 spectrum at scales smaller than a whereas it can be easily shown that 
the corresponding IK arguments lead to a k~ 5 ^ 2 law. 
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4 Karman-Howarth theorem for LAMHD—a 
in 2D and 3D Elsasser variables 

In 1938, Karman and Howarth [21] introduced the invariant theory of isotropic 
hydrodynamic turbulence, and derived from the Navier-Stokes equations an 
exact law relating the time derivative of the two-point velocity correlation 
with the divergence of the third-order correlation function. Later, this result 
was generalized to the MHD case by Chandrasekar [213, and recently written 
in terms of Elsasser variables [22] • For LANS— a in the fluid case it was de- 
rived in J2H] • The relevance of the Karman-Howarth theorem for the study of 
turbulence cannot be underestimated. As a corollary, rigorous scaling laws 
in the inertial range can be deduced. In this section we will generalize these 
results to the LAMHD—a case. 

For the sake of simplicity, we will consider the case r\ — v — 0, the 
dissipative terms can be added at any point in the derivation. Also, we will 
use a = au- We start writing the LAMHD—a equations using the Elsasser 
variables 

z ± = u±B, z s ± = u s ±B s . (23) 
Applying the Helmholtz operator to eq. (j2J), we obtain 

(1 - a 2 A) {d t B s + u s • VB S - B s ■ Vu s ) = 0. (24) 

Now we add and subtract eqs. (fTTj) and (j21j). Using eqs. (J22|) we obtain 
equations for the evolution of z s ± , 

(1 - a 2 A) (dtz* + z s T ■ Vz/ + W s ) = -« 2 div r, (25) 

where the stress tensor divergence div r in terms of the Elsasser variables is 

div r = - div (Vz s + ■ Vz s ~ + Vz s + ■ Vz s ~ T - Vz s +T • Vz s " 

+Vzr ■ Vz s + + Vz s ~ ■ Vz s +T - Vz s - T ■ Vz s + ) . (26) 
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This stress divergence may be rewritten equivalently, as 

divr = 2div(A + -£- + A~ - £+) + Vn 

with S ± = i(vz ± + (Vz ± f), A ± = ^(vz ± -(Vz ± ) T ), 

with the same additional pressure II as before. 

We could repeat all the derivation to obtain an equation for the evolution 
of from eqs. (JTJ) and (j2U). Instead, starting from eq. (|23|) . using eqs. (jlj) 
and (|TH) we obtain 



where 



9 t z ± + z s T ■ Vz ± + W = a 2 div q ± (27) 



? ± = _ ( Vz/ ■ Vz s T + Vz/ • Vz s tT + Vz s ±T • Vz A T 



2 

-Vz s T • Vz/ - Vz s T ■ Vz s ±T + Vz s tT • Vz/) . (28) 



Note that equations (j2*3j) and (|2*7j) make explicit the fact that Alfven 
waves u = ±B, u s = ±B S are exact nonlinear solutions of the LAMHD—a 
equations. For an Alfven wave either z + or z (as well as the corresponding 
field z^) is zero. In this case, all nonlinear terms are zero and verification 
of the solution follows. 

In Cartesian coordinates, we can write equations (|25|) and (|27|) in com- 
ponents 

dtzf + d k (zfz s ^ k + tW* - a\ ±k ^ = (29) 
d t z' s f + d' k {zfz^ k + + a 2 g a * r>fj = 0, (30) 



the prime denotes that the variables are evaluated at x', and 

-r/a 

is the Yukawa potential. The Green function of the Helmholtz operator is 



9c = T—r (31) 
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given by 

9a*r* = J ^(|x'-x|)rf(x')dV, (32) 
and the components of the stress tensors r and q ± are 

T* = ~ (d 3 Z s td k Z- ] + d jZs t&Z- k - d iZ +d k Z s ~ j 

+d jZs -d k z s +3 + d 3 z s -&z+ k - d lZs -d k z s + ^ , (33) 

= l(d j z a t&z* + d j z a t#z* h + d i z a f&z a : V 

-djz a fd k z s ±j - d 3 z s ?cPz a ±k + d i z a Jd k z a ±j ) . (34) 



Multiplying eq. by z'^, eq. (|3U|) by zf, and adding the result yields 

d t (z±z'f) = ^/( z ^ z ^- a ^ k ) z 'f\ + l T (nz'f5 k -n'z^5 k 



1 \ Z i Z sj J Q r k \ \ Z i Zs a *> ' » J ~sj I ' Q rk y~sj "i "s~i "j 



— ( ( ?' ± ~' 

dr k \\ s i 



z n - z„ + a 



9a*r';)z^), (35) 



where we used homogeneity 



dr k ^ ^ dx' k ^ ^ dx k ^ ^ 

Now, we can make the equation symmetric in the indices i,j adding the 
equation for d t (zi z' s f ). We use homogeneity 

QiQ aj Q' s k + W'sil's*) = ~ (l'i<ls j <ls k + q'^qJ") , (37) 
and define the tensors 

QS = U<t + zt<t), (38) 



r% = {(ztz' s t + zfz' s t + zfz s f + zfz s t)z a ^), (39) 

= ((K4-*Zt) 5 ! + (tt-*Zt) 5 i)> ( 4 °) 



> y = ( r i^ + T "i 2: s i +9a*T aj zf + g a *T ai zf 
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We can drop n^- because the terms with the pressures tt and n' a vanish 
everywhere, as follows from the usual arguments of isotropy [21] . Finally we 
obtain 

This is the LAMHD—a version of eq. (3.8) in [21]. In the case a = this 
equation is also a linear combination of eq. (43), (50), and (56) in [22]. More 
Karman-Howarth equations can be written for different combinations of z ± 
and zf. 

Since Qfj and T^\j are symmetric and divergence free in their indices i and 
j, <S ± ij must be symmetric and divergence free in i and j. But the Elsasser 
variables are combinations of vectors and pseudovectors. Therefore, Q ± 
is a combination of tensors and pseudotensors. We can define a tensor as 

Q=Q ± + Q T , (43) 

and a pseudotensor as Q — Q T . We will continue using only the tensor Q, 
the results can also be obtained for the pseudotensors using the expressions 
in [21]. We also define T = T± + T T and S = S ± + S T . 

Imposing isotropy and from incompressibility, Q can be written as [2Z] 

Qij = cmliQne^) = -(d+ l)Q8ij + rQ' - , (44) 

where the curl is taken with respect to the third index (j), Q = Q(r,t) is 
a scalar function, e is the Levi-Civita pseudotensor, and d is the number of 
dimensions. Here, Q' = d r Q. 
In the same way, we can write 

7~J = curl [T{ne jk ir l +r j e m r 1 )] 

= -T'r irj r k - (rT' + dT) (r^J + r$ ) + 2T5 ijr k . (45) 

The tensor S takes the same form with scalar function S(r, t). Note that S is 
the isotropic sub-a-scale stress tensor in the LES formulation of LAMHD—a. 
Now we compute the divergence of these tensors. In three dimensions 

A-T k = cnTl[(rr + 5T)e tjl r l ], (46) 
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and the divergence for S takes the same form. Replacing eqs. (|44j) and (|46|) 
into eq. (|42j) we finally obtain 

Theorem 1 (Karman-Howarth Theorem for LAMHD—a ) 
The exact LAMHD—a model relation \^ 6 4) for homogeneous isotropic statis- 
tics implies the isotropic tensor relation in three dimensions 

d 4=(4 + s)(T 2 -^)> (47) 



(T 2 - a 2 S) . (48) 



dt \ dr 
and in d dimensions the general result is 



dQ 
dt 



r|- + (d + 2) 



This is the generalization of the Karman-Howarth equation for LAMHD—a 
(two more equations can be written for different combinations of the tensors 
and pseudotensors) , without the dissipation. When B = this equation is 
also eq. (3.16) in [2E]. When a = 0, this is equivalent to the Karman-Howarth 
equation for the Elsasser variables as derived in [22], or a combination of 
equations (49) and (53) in [25] . 

Therefore, all equations in (221 follow for a/r <C 1. This result confirms 
that the alpha-model preserves the properties of MHD for separations larger 
than r ~ a. For r > a, the scaling of structure functions and the relation 
between second and third order functions hold. 



Corollary 2 (Kolmogorov Theorem for LAMHD—a ) 

Introducing the flux d t Q = —2e a /d with e a = e+ + e~ (the energy injection 
rate for each Elsasser variable) in eq. J^&| ) and integrating in the inertial 
range yields 

-5(5T2) £ « = ( T - Q ' 2s ) • < 49) 

where T and S are defined in equations \39^l and \41\j . 

Note a multiplicative factor compared with the usual expression from Kol- 
mogorov, related to the relation between autocorrelation functions and struc- 
ture functions in isotropic turbulence in d dimensions. For a/r < 1 this 
equation reduces again to the MHD results. Note also that structure and 
autocorrelation functions in LAMHD—a involve one unsmoothed field and 
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one smoothed field if quantities are of second order, and two smoothed fields 
if quantities are of third order. In the following sections, we will use this con- 
vention and all structure functions for LAMHD—a will be written as they 
follow from the expressions of the tensors Q and T. 

5 Energy dissipation anomaly 

The Karman-Howarth (KH) theorem for fluid turbulence [21] gives the ex- 
act analytical relation between the time rate of change of the second-order 
two-point velocity correlation function and the gradient of the third-order 
two-point velocity correlation function, as derived from the Navier-Stokes 
equation for homogeneous, isotropic turbulence. 

Kolmogorov JH] used the structure function form of the KH equation, 
to show - for homogeneous, isotropic and stationary turbulence, in the limit 
v — > of vanishing kinematic viscosity - that the Navier-Stokes equations 
lead to an exact relationship between the third-order structure function and 
the energy dissipation rate, e, which scales linearly in the separation, r. This 
is Kolmogorov's famous "four-fifths law." 

By assuming self-similarity of scales in the inertial range Kolmogorov 
then was able to deduce, in steps that essentially amount to dimensional 
analysis, that the second-order structure function must scale as r 2//3 and that 
consequently the energy spectrum (which is essentially the Fourier transform 
of the second-order structure function) must scale as k~ 5 ^ 3 . As noted in |2~H] . 
Kolmogorov's four-fifths law 

is one of the most important results in fully developed turbu- 
lence because it is both exact and nontrivial. It thus constitutes 
a kind of 'boundary condition' on theories of turbulence: such 
theories, to be acceptable, must either satisfy the four-fifths law, 
or explicitly violate the assumptions made in deriving it. 

The two key assumptions in Kolmogorov's derivation of the four-fifths law 
are that (1) an inertial range exists in which the flow is self-similar and (2) 
the energy dissipation rate does not change as one takes the limit v — > 0. 

The equivalent of the KH equation was derived for the LANS— a model 
in Since the model relates the Helmholtz smoothed velocity u s , to the 
unsmoothed velocity u, the appropriate correlation functions that emerge 
are the second- and third-order two-point correlation between u s and u. 
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Upon following Kolmogorov's analysis for isotropic inertial range statistics, 
the corollary to the LANS— a KH-equation is that solutions of the LANS— a 
equations possess two regimes of scaling, depending on whether the separa- 
tion distance r is greater, or less than the size a. For r > a, the third-order 
correlation scales like r, thereby recovering Navier-Stokes behavior. In con- 
trast, for r < a the third-order correlation scales like r 3 . If self-similarity 
is then assumed one finds for r > a that the second-order correlation scales 
like r 2//3 , again recovering Navier-Stokes behavior. However, for r < a the 
second-order correlation scales like r 2 . Correspondingly, the power spectrum 
for the smoothed velocity u has two regimes, with a transition from k~ 5 ^ 3 
for ka < 1 to k~ 3 for ka > 1. Thus, the KH-theorem for the LANS— a 
model derived in is consistent with the spectral scaling results found for 
it in [22] by dimensional arguments. We shall apply similar reasoning to the 
LAMHD—a model in two dimensions. 

Differences from MHD turbulence theory for r < a 

The second term in the il —2/d(d+2) Law" in equation (J49)) (the a 2 S term on 
the right side) is reminiscent of the quantity that appears in the correspond- 
ing "—2 Law" for enstrophy cascade in 2D turbulence. The latter expression 
contains two powers of enstrophy and one power of velocity. For example, 
see the Appendix B of where this identity for 2D turbulence is derived 
in detail. 

Likewise, the a 2 S term in ([4"9"|) for LAMHD—a in 2D contains two pow- 
ers of gradients Vz s and one power of z s without gradient. Consequently, 
this should be the dominant term (compared to the first T— term) for small 
separations, when r < a. If in addition the LAMHD—a flow is self-similar, 
the dominance of the a 2 S term in (J49|) when r < a implies a scaling re- 
lation for the second-order structure functions. Following [T^j as amplified 
by [2E]> let the longitudinal difference 5z(x, r) satisfy the scaling relation 
5z(x, Ar) = A h <5z(x, r) for all x and for increments r and Ar small compared 
to a. By dimensional analysis, [<S(Ar)] = [(5z) 3 /r 3 } = [«S(r)]. Consequently, 
3h — 3 = and h = 1 for small scales r < a in a self-similar LAMHD—a 
flow. This means the second-order structure functions follow r 2 scaling for 
r < a in such a flow. This r 2 scaling implies a k~ 3 law for the spectral 
density of smoothed kinetic and magnetic energy in that range for the 2D 
LAMHD—a model. Thus, one finds a self-similar k~ 3 "enstrophy-like" cas- 
cade, in agreement with similar considerations of (22] for the kinetic energy 
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spectral density in the 3D LANS— a model. 

Implications of the A;" 5 / 3 — > k~ 3 spectral scaling transition for the 
LAMHD—a dissipation anomaly. The second term in the u —2/d(d + 2) 
Law" in equation (|%9[) modifies Kolmogorov's four-fifths law at small sepa- 
rations (r < a), provided one may assume constancy of total LAMHD—a 
energy dissipation as v — > 0. This is the "energy dissipation anomaly" for 
the LAMHD—a model. A technical argument using embedding theorems for 
Besov spaces first introduced in [3U] implies that constancy of total energy 
dissipation may hold as v — > for a turbulent fluid in two dimensions, pro- 
vided its L 2 power spectrum is not steeper than k~ 4 . The k~ 3 spectrum for 
ka > 1 is not too steep; so the roll-off k~ 5 ^ 3 — > k~ 3 in the LAMHD—a power 
spectrum is consistent with the necessary condition for possessing such an en- 
ergy dissipation anomaly. Hence, the k~ 3 behavior in the L 2 power spectrum 
of the LAMHD—a model for ka > 1 and the corresponding modification for 
separations r < a of Kolmogorov's four-fifths law derived in are both 
consistent with the assumption of constant dissipation of total kinetic energy 
as the Reynolds number tends to infinity. 

Implications of the spectral scaling transition for IK scaling instead 
of Kolmogorov scaling. While in hydrodynamic turbulence the spectral 
transfer is believed to be a local process, in MHD turbulence Iroshnikov and 
Kraichnan proposed that the spectral transfer is governed by nonlocal Alfven 
wave interactions (see [21] for a study of non-locality of transfer in MHD). 
If a large scale magnetic field is present, then this field acts as a guide field 
to the fluctuations, turning them into Alfven waves. Kraichnan proposed 
that Alfvenic propagation limits the nonlinear interaction responsible for 
the transfer of energy to smaller eddies in the absence of magnetic fields. 
As a result, the IK spectrum is shallower than K41. In LAMHD—a , the 
k~ 3 spectrum at scales smaller than a turns into a shallower k~ 5 ^ 2 when 
the IK hypotheses are used. Note that the weakening of local interactions 
due to Alfven waves holds in LAMHD—a at the smoothed scales, and the 
energy dissipation anomaly should also be captured. This is also in agreement 
with equations and ([270. which show that Alfven waves are also exact 
nonlinear solutions of the LAMHD—a equations. 
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6 Numerical results for intermittency and scal- 
ing anomaly 

In this section, we compare intermittency in LAMHD—a to that of direct 
numerical simulations (DNS) of the MHD equations, regarded as true at a 
given Reynolds number. Intermittency is associated both with the presence 
of strong localized structures and with the existence of strong non-Gaussian 
wings in the probability distribution functions. We have previously inves- 
tigated the latter in ^1] and we concentrate here on the strong localized 
structures giving rise to deviations from universality, as can be studied by 
examining high order statistical moments, such as the structure functions. 

We define the longitudinal structure function of a field f as = 
(|^/l| p ) where 8fi = (f(x + l) — f(x)) ■ 1/7 is the longitudinal increment of f. 
In the inertial range between the large energy- containing scales and the small 
dissipative scales, the structure functions are assumed to vary in a self-similar 
manner, ~ fa ■ As previously mentioned, in isotropic and homogeneous 

turbulence the structure functions can be related to the correlation functions 
discussed in Section 0] K41 phenomenology predicts Q = p/3, while IK 
gives = p/A. The anomalous departure of the exponents ( p from these 
linear scaling laws is a measure of intermittency-induced deviations from 
universality. 

To numerically solve the MHD and LAMHD—a equations we will use a 
parallel pseudospectral code as described in [TT]. In two dimensions, the ve- 
locity and magnetic field are expressed as the curl of a scalar stream function 
\P and a one component vector potential A z , respectively: 



where \1> = (1 — a 2 V 2 )\l/ s , and A z — (1 — a 2 M W 2 )A Sz . In terms of these 
quantities, the 2D MHD equations may be expressed as 



u = V x (*z), u s = V x (\M 



(50) 



B = V x (A z z), B s = V x (A Sz z) 



(51) 



d t V 2 ^ = V 2 ^] - [A t , V 2 A Z ] + zA7 2 V 2 # 



(52) 



d t A z = [$,A z }+r)V 



2 A 



Z1 



(53) 



18 



where 



[F, G] = d x Fd y G - d x Gd y F 



(54) 



is the standard Poisson bracket. The LAMHD—a equations (012)) modify 
this two dimensional structure by introducing smoothed variables as 



In the following subsections we test the LAMHD—a model against MHD 
results (for which a = au = 0) for freely decaying turbulence with the 
same initial conditions, dissipation and time-stepping, and also for forced 
turbulence where we have averaged statistics over 189 turnover times taken 
from 9 experiments with distinct seeds for the random forcing, resulting in a 
data set of ~ 2 • 10 s points. 

6.1 Forced simulations 

In this subsection we consider forced turbulence with rj — v — 1.6 x 10 -4 . 
Four sets of simulations were carried out, one set of MHD fully-resolved 
simulations with 1024 2 grid points, and three sets of LAMHD—a simulations, 
with 512 2 grid points and a = au = 6/512, with 256 2 grid points and 
a = au = 6/256, and with 256 2 grid points and a = au = 6/128. Note that 
the 256 2 LAMHD—a simulation with a = au = 6/128 could be carried out 
with a 128 2 spatial resolution (see e.g. criteria for the selection of the value 
of a and the linear resolution in ^1]). The reason to keep the resolution fixed 
at 256 2 is to preserve the amount of spatial statistics, crucial to measure high 
order exponents as will be shown later. 

Both the momentum and the vector potential equations were forced. The 
expressions of the external forces were loaded in the Fourier ring between 
k = 1 and k = 2, and the phases were changed randomly with a correlation 
time At = 5 x 10 -2 . Averaged over space, the amplitudes of the external 
forces were held constant to Fm = 0.2 in the vector potential equation, 
and Fk = 0.45 in the momentum equation. The systems were evolved in 
time until reaching a turbulent steady state, and then the simulations were 
extended for 21 turnover times. Over this time span, 21 snapshots of the 
fields from each run were used to compute the longitudinal increments. As 



«9 t V 2 # = [tf a , V 2 ^] - [A 8z , V 2 A Z ] + i/V 2 V 2 * 



(55) 



d t A Sz = [q s ,A Sz ]+ v V 2 A z . 



(56) 
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previously mentioned, each set of simulations (DNS, and LAMHD—a with 
different spatial resolutions) comprises nine runs with the same viscosity and 
diffusivity but different series of random phases in the external forcing, to 
have enough statistics to determine the scaling exponents up to eighth order 
with small error bars. The total number of points was ~ 2 ■ 10 8 for the DNS 
and ~ 1.2 • 10 7 points for the 256 2 LAMHD—a simulations. 

During these intervals, the integral Reynolds number based on the r.m.s. 
velocity fluctuates around 2200. The normalized correlation coefficient be- 
tween the velocity and the magnetic field is 20% with a standard devi- 
ation of 16% within the set of nine runs, and its unsigned counterpart 
2(|u ■ B|)/(|u| 2 + |B| 2 ), is ~ 29% ± 12%. The ratio of the integral scale 
to the Taylor scale computed on the r.m.s. fields fluctuates around 10 for 
all the simulations. The ratio of magnetic to kinetic energies is ~ 2 for all 
runs. Finally, the Kolmogorov dissipation wavenumbers k v = ((uj 2 ) /V 2 ) 1 / 4 
and krj = ((j 2 ) /7/ 2 ) 1 / 4 fluctuate around 330, values substantially larger than 
the largest resolved wavenumbers k a ~ 1/a in all LAMHD—a simulations, 
by virtue of the model. 

Average omni-directional spectra for magnetic and kinetic energies over 
these 189 turnover times are shown in Figure H All spectra display an iner- 
tial range, and the LAMHD—a simulations are able to capture the spectral 
behavior up to k « 1/a. For k > 1/a, theoretical arguments suggest a k~ 3 
spectrum for the alpha model. To observe this spectrum, however, large 
values of a would be required. This is inconsistent with the use of the al- 
pha model as a subgrid turbulent model, and lies beyond the interest of the 
present work. 

K41 theory predicts 

e£(Q (57) 

which follows from the assumption that the statistical properties of the field 
are self-similar in the inertial range, which will be identified here as the scales 
for which the relation 

C3 = 1 (58) 

holds. The existence of scaling (J57)) has been extensively verified for the 
hydrodynamic case. Starting from the assumption of self-similarity in the 
inertial range, we can then postulate the validity of (|57|) at arbitrarily high 
order, p. From an experimental standpoint, the amount of data used in 
calculating &p(l) determines the highest order for which we can observe 
a scaling in the inertial range. Above this order, the assumption of self- 
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Figure 1: Spectra averaged over 189 turnover times. 1024 2 MHD is the solid 
line, 512 2 LAMHD—a is the dotted line, 256 2 LAMHD—a is the dashed line, 
256 2 LAMHD—a with a = «m = 6/128 is the long-dashed line (hereafter 
indicated in figures as '256 2 2a'), k~ 5/3 (K41) is the dash-dotted line, and 
/j-3/2 (j^q j g dash-triple-dotted line; the K41 and IK slopes are shown 
for reference. The vertical lines indicate the wavenumbers corresponding to 
the lengths a for all LAMHD—a simulations. Panel (a) is magnetic energy, 
Em, versus wavenumber k, and panel (b) is kinetic energy, E K vs. k. 
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similarity allows us to fit a scaling law to our data in any event. But in 
the absence of sufficient statistics, such a fit can be poor and the error bars 
rather large (see, e.g., FigureEJand Figure El (a) to be discussed shortly). 
The Extended Self-Similarity (ESS) hypothesis [32] , proposes the scaling 

6£(Q ~ , (59) 

which is found to apply to a much wider scaling range than the inertial 
range. Here, the scaling range is determined by the observed scaling for low 
orders and a chief benefit is increased statistics to compute more accurate 
exponents at higher order. For the case of MHD, Ref. [22] proposes to replace 
63(f) with the third-order, mixed structure functions, L ± {1). As in the ESS 
hypothesis, these scales L ± (l) may provide better independent variables (as 
opposed to length I) against which to determine the scaling exponents for 
MHD. From the KH theorem for an incompressible, non-helical MHD flow, 
they find 

(5ztmz ± (\)\ 2 ) = ~e ± l, (60) 

where |fc ± | 2 = (5z±) 2 + (Sz$) 2 , 5z$ are the transverse increments, d is the 
space dimension, e + and e~ are the energy dissipation rates for \{z + ) 2 and 
I (z~) 2 respectively, and angle brackets indicate as usual spatial averages [22] • 
These results for the third-order structure functions are exact and can 
be used to compute more accurate anomalous scaling exponents of structure 
functions of higher order. Due to cancellation problems (linked with having 
limited statistics), absolute values are often employed; we also find linear 
scaling in this case, viz.: 

L ± (l) = (I^H^I 2 ) ocZ . (61) 

As follows from the expressions given in Sec. [I] and the invariants found 
for both MHD and LAMHD—a [TTl I14j . when making comparisons between 
DNS and model runs, we substitute the H\ norm, (|H| 2 ) = (|u-u s |) [TBI 133] . 
for the regular L 2 norm, (\u\ 2 ) = (|u • u|), whenever we consider quantities 
for the LAMHD—a model. The Karman-Howarth theorem for LAMHD—a 
is essential to this study of intermittency in that it allows us to define the 
structure functions for LAMHD—a; it also identifies the flux relation that 
scales linearly with / for application in MHD of the ESS hypothesis. Accord- 
ingly we determine the relative scaling exponents, by using Eq. (|8T|) for 
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the third-order, mixed structure function, L + {1) = {\Sz L \ \5z + \ 2 ) for MHD 
and L+(l) = (\5z; L \\\5z + \\l) for LAMHD-a, 

& f P (l) ~ HM 4 ■ ( 62 ) 

Figure 121 shows the third-order mixed structure function L^(l)/l as a 
function of I. We find, contrary to what is reported in J7|, that the relation 
(16 lj) has a identifiable range of validity, as can be seen in the figure by 
comparison with the solid straight line denoting the computed slope L + (l) ~ 
This range of validity is identified as the inertial range, 27r/20 < I < 
27r/10, indicated by dashed vertical lines. The LAMHD—a runs display the 
same scaling as the MHD simulation, and departures are pronounced only 
for scales approaching and smaller than a (for the 256 2 runs, a ~ 0.15, 0.29 
and for the 512 2 run, a f=s 0.07 as indicated by dotted vertical lines). Note 
that the results have been scaled by the mean value of LtM). As the average 
energies of the runs are disparate, this improves the ease of comparison. The 
same behavior is observed for L7M). 

The scaling of the third-order structure function Sjj~ versus for the 
Elsasser variable z + is shown in Figure El (a) as well as a compensated plot 
versus / in Figure|3](b). We see very little contamination at scales larger than 
a. This contrasts with hyperviscosity which is known to cause an enhanced 
bottleneck in Navier-Stokes turbulence, which corrupts the scaling of the 
larger scales for structure functions of order two and higher The MHD 
case has not been studied in this context, and the presence of an inverse 
cascade of magnetic helicity might also exacerbate the problem. 

The LAMHD—a simulations show similar scaling to what is found in 
DNS of MHD. A solid line indicates the best fit to the DNS data, & = 
1.01 ± 0.08 using the ESS hypothesis (the ESS scaling range is indicated 
by dotted vertical lines); note that all errors presented and shown in the 
figures correspond to 3a where a is the standard deviation. The much smaller 
inertial range is indicated by dashed vertical lines and arrows indicate lengths 
a. This is the main benefit of using the ESS hypothesis. While in the 1024 2 
simulation there is enough statistics to measure the scaling exponents £ p in 
the inertial range, when we use the LAMHD—a equations to reduce the 
computational cost, the amount of spatial statistics is drastically reduced 
(e.g. by a factor of 16 in the 256 2 runs). The ESS hypothesis allows us 
to extend the range where the £ p exponents are computed, giving a better 
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Figure 2: Third-order, mixed structure function, L^(/)// 0,99 , versus 1, for 
forced runs of turbulence averaged over 189 turnover times. Results are 
scaled by the mean value of LfM) for easier comparison. Labels are as in 
Fig. [TJ The best fit to the MHD data, L + ~ /° ", is indicated by the solid 
straight line. The inertial range where this fit is made is indicated by dashed 
vertical lines and dotted vertical lines indicate the lengths a for the 512 2 and 
256 2 simulations. 
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Figure 3: Third-order structure functions for z + : Panel (a) versus Lt^ 

and panel (b) versus /, computed over 189 turnover times. Labels 

are as in Fig. ^ The solid line corresponds to the best fit to the DNS 
data, &3 = (L + ) 101 (the dash-dotted lines represent the 3a error). The ESS 
hypothesis range where this fit is made is indicated by dotted vertical lines 
and dashed vertical lines indicate the inertial range. Arrows indicate the 
several lengths a used in the simulations. 




Figure 4: Eighth-order structure functions for z + : Panel (a) (Sg versus 

and panel (b) /l^ versus /, computed over 189 turnover times. Labels 
are as in Fig. The solid line corresponds to the best fit to the DNS data, 
= (L + ) 1,63 . The ESS hypothesis range where this fit is made is indicated 
by dotted vertical lines and dashed vertical lines indicate the inertial range. 
Arrows indicate lengths a. 
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estimation and smaller error bars. As an example, in Fig. 0] (a) we show 
the scaling of the eighth-order structure function (3g versus L~f\ as well as 
a compensated plot versus I in Figure |U (b). The ranges corresponding to 
the inertial range and ESS are also indicated. For the 256 2 runs, we cannot 
observe a scaling at this order. The error of a scaling computed from the 
assumption of self-similarity is excessively large (see Figure El (a)). From the 
ESS hypothesis a better estimation of the (postulated) scaling at order eight 
can be made (see Figure El (b)). 

Figure El compares the scaling exponents, for the DNS runs and the 
three sets of LAMHD—a runs. Figure El (a) is for exponents computed only 
over the inertial range. Notice that the LAMHD—a runs capture the scaling 
for the low-order moments (p < 4). For higher-order moments (beginning 
at p = 5), the drop in the scaling exponents for the 256 2 results (with a = 
atM = 6/256) and the large error bars are indicative of insufficient statistics. 
The advantages of ESS are clearly seen by comparison with Figure El (b), 
which shows the scaling exponents for all sets of runs employing the ESS 
hypothesis. In both figures, the She-Leveque (SL) formula [33] modified for 
the MHD case [21] is shown as a reference, 



From these results, we conclude that LAMHD—a captures the intermittency 
of the DNS runs up to and including the eighth-order moment (to within the 
errors of our statistics). The size of the error bars make it difficult to draw 
further conclusions from the data. 

The anomalous scaling results for the DNS runs are shown in Table [0 
Though our goal here is to test LAMHD—a against DNS, we remark briefly 
on the correspondence between our scaling exponents and other studies. In 
opposition to the findings of [2Ej, we find £3 ~ 1 but £f > 1. We note, 
however, that the forcing in was tailored to maintain at a constant level 
all Fourier modes with k — 1 while our forcing is random with a constant 
amplitude between k — 1 and k = 2. As can be seen from Table ^ our 
results for are in good agreement with jTJj, for decaying turbulence, and 
with j37j, for forced turbulence. As [T7j suggests, the scaling exponents, as 
inertial range properties, may depend on the character of the driving due to 
non-local processes in the cascade dynamics connected with the Alfven effect 
(see also anc ^ as ^ e f- Ell points out, such an analysis can be sensitive 




(63) 
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Figure 5: Structure function scaling exponent: versus p, computed over 
189 turnover times. 1024 2 MHD are the pluses, 512 2 LAMHD-a are the 
diamonds, 256 2 LAMHD—a (a = 6/256) are triangles, and 256 2 LAMED— a 
(a = 6/128) are the squares. The dashed line indicates K41 scaling, dotted 
line indicates IK scaling, and the solid line is the prediction using the modified 
She-Leveque formula (see text). Panel (a) is computed over the inertial 
range. Panel (b) is computed utilizing the ESS hypothesis. 
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Table 1: Relative scaling exponents (together with 3a errors in computing 
the slope in parenthesis) computed from the 9 MHD runs over the inertial 
range, £, and utilizing the ESS hypothesis, £,ess- 



V 




(ess 


(ess 


Z.ESS 


Z.ESS 


1 


.43(10) 


.43(05) 


.42(04) 


.37(06) 


.43(04) 


2 


.75(13) 


.76(06) 


.76(04) 


.67(08) 


.76(04) 


3 


.99(14) 


1.01(08) 


1.01(03) 


.92(08) 


1.01(04) 


4 


1.16(14) 


1.19(09) 


1.20(04) 


1.10(08) 


1.16(02) 


5 


1.29(11) 


1.34(10) 


1.34(07) 


1.25(08) 


1.27(04) 


6 


1.37(09) 


1.45(11) 


1.45(09) 


1.36(05) 


1.33(06) 


7 


1.43(13) 


1.55(13) 


1.53(12) 


1.44(11) 


1.38(07) 


8 


1.46(22) 


1.63(15) 


1.60(14) 


1.50(13) 


1.43(08) 



to several parameters such as the ratio of kinetic to magnetic energy or the 
amount of correlation between the velocity and the magnetic fields. 

Figure El shows the scaling exponents for the velocity and magnetic fields, 
as well as for the other Elsasser variable z~. The anomalous scaling is again 
matched by LAMHD—a up to and including eighth-order. The results from 
the MHD and LAMHD—a simulations are also in good agreement with Ref. 
|37j . Note that the magnetic field is more intermittent than the velocity 
field (in the sense that the scaling exponents deviate more from a straight 
line), as previously found in numerical simulations it may be related to 
the fact that in MHD, nonlinear interactions are more non-local (in Fourier 
space) than for fluids jUJ. This well known feature of MHD turbulence is 
also properly captured by the LAMHD—a equations. The average (over all 
fields) of the 3a errors of the eighth order scaling exponent is 0.13 for the 
DNS on a 1024 2 grid; it is 0.15 and 0.16 for LAMHD-a on a 512 2 and 256 2 
grid, respectively. To further test the convergence of our statistics, we re- 
duced the amount of data used to compute the scaling exponents for the 
512 2 runs by a factor of 4 (which gives the same amount of statistics than 
the 256 3 LAMHD runs) and determined an average 3a error of 0.17. While 
these results confirm our convergence as the amount of statistics is increased, 
they also highlight the rather low decrease in error with increased compu- 
tational effort. Accordingly, the computational burden for more accurate 
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determination of high-order statistics prohibits further analysis of the data. 
6.2 Decaying simulations 

In this subsection we briefly discuss simulations of free decaying turbulence 
using both the MHD and the LAMHD—a equations. The results are similar 
to the ones presented in the previous subsection for forced turbulence. How- 
ever, since no turbulent steady state can be defined in freely decaying runs, 
the amount of statistics is reduced as only a few snapshots of the velocity and 
magnetic field during the peak of mechanic and magnetic dissipation can be 
used to compute structure functions. To partially overcome this problem, we 
will discuss simulations with higher spatial resolution than the ones presented 
in the forced case. 

A fully resolved 2048 2 MHD run was made using v = t] = 10~ 4 , as well as 
a 1024 2 LAMHD—a run with a = a M = 6/1024 and a 512 2 LAMHD—a run 
with a = a^f = 6/512. The initial velocity and magnetic fields were loaded 
with random phases into the rings from k = 1 to k = 3 in Fourier space. The 
initial r.m.s. values of u and B are equal to unity. No external forces are 
applied and the system decays freely as a result of the dissipation. Under 
these conditions, the LAMHD—a equations have been shown to reproduce 
the time evolution of the magnetic and kinetic energy, as well as the evolution 
of the spectra and other statistical quantities [TT] . 

The magnetic and kinetic energy spectra between t = 3 and t = 6 in 
units of the eddy turnover time (the time for which a quasi-steady state 
is observed) display an inertial range with an extent of approximately one 
decade in Fourier space, from k ~ 3 up to k ~ 30. As a result, one snapshot of 
the fields in this range in time was used to compute the longitudinal structure 
functions. The Kolmogorov kinetic and magnetic dissipation wavenumbers 
k v and k v peaked at a value of 470, larger than the filtering wavenumber k a ~ 
1/a in all the LAMHD—a simulations. We note that for large wavelength 
component behavior up to k ~ k a , both LAMHD—a simulations accurately 
reproduced the omni-directional spectra for the magnetic and kinetic energies 
as was shown previously for the forced runs. 

Both LAMHD—a runs preserve the scaling of the longitudinal structure 
function exponents observed in the MHD simulation. As an example, Fig. 
13 shows the £ p exponents for the z + Elsasser variable for the MHD and 
LAMHD—a simulations using the ESS hypothesis. Note that the £ p expo- 
nents of the three simulations lie within the error bars, and the three simula- 
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Figure 7: Scaling exponent as a function of the order p, for the z + Elsasser 
variable in simulations of freely-decaying turbulence. The MHD simulation, 
the 1024 2 LAMHD—a run, and the 512 2 LAMHD-a run are denoted by 
crosses, diamonds, and triangles respectively. The K41 scaling is the dashed 
line, the IK prediction is the dotted line, and the solid line corresponds to 
eq. (pj) . 

tions show departures from the self-similar K41 or IK scaling. For values of p 
larger than 6, effects associated with the limited amount of spatial statistics 
can be observed in all the runs. 



7 Discussion and outlook 

Sufficient resolution for studying directly high Reynolds number flows as 
encountered in geophysics and astrophysics is today well beyond technological 
limits. Closures such as the Lagrangian-averaged alpha-model can reduce the 
computational burden by reducing the resolution requirements. However, to 
be used as a model of hydrodynamic or magnetohydrodynamic turbulence, 
or for applications in astrophysics and geophysics, detailed knowledge of the 
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ability of the LANS— a or LAMHD—a equations to capture key features of 
turbulent flows is required. 

The LANS— a and LAMHD—a equations have been tested against direct 
numerical simulations in a variety of problems (see e.g. jHH] for neutral fluid 
studies and [TT1[T21[I3] for studies in conducting fluids). Most of these works 
compared the time evolution of ideal invariants for forced and free decaying 
turbulence, as well as the evolution of energy spectra. Also, some statistical 
comparisons were performed (e.g. studying the behavior of probability den- 
sity functions). In this work, we apply a more stringent test to these models. 
Intermittency is a well known feature of turbulent flows, associated with the 
existence of strong events localized both in space and time. Intermittency 
can trigger large scale events, affect the transport coefficients, or give rise to 
corrections in the turbulent scaling. As a result, whether a model can cap- 
ture the statistics of intermittent events is of utmost importance to model 
astrophysical or geophysical flows. The study of intermittency also requires 
computation of high order statistics, thereby extending previous comparisons 
between DNS and a-models. 

An extension of the Karman-Howarth theorem (KH-a) was proven for the 
LANS— a equations in Ref. [211] • As a corollary of this theorem, Kolmogorov 
four-fifths law and Kolmogorov's energy spectrum can be derived for the 
LANS— a equations at scales larger than a. In this work, we extended the 
KH-a theorem to the LAMHD—a case, proving as a result that for scales 
larger than a the LAMHD—a equations satisfy the scaling laws previously 
known for MHD turbulence [22 123 • This is an important result, since MHD 
turbulence involves two coupled fields (the velocity and magnetic fields) and 
can display different power laws in the inertial range according to the regime 
of interest. While LES often impose a particular regime and a power law 
El M HO], the LAMHD—a equations are shown to satisfy the general 
scalings satisfied by the MHD equations without any hypothesis about the 
scaling followed in the inertial range. 

The extension of the KH-a theorem to the LAMHD—a case also allows 
us to define correlation and structure functions in the a-model. With these 
functions, the analysis of anomalous scaling and intermittency can be per- 
formed. Numerical simulations were carried out both for freely decaying and 
for forced two dimensional MHD turbulence, solving directly the MHD equa- 
tions, and employing the LAMHD—a equations at 1/2 and 1/4 resolution (a 
case equivalent to 1/8 of the DNS resolution was also considered for forced 
turbulence). In the forced runs, we have averaged statistics over 189 turn- 
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over times (and up to ~ 2 • 10 8 points) to test if the LAMHD—a equations 
reproduces intermittent turbulent behavior. The scaling of the third-order 
structure function was tested and linear scaling with length (down to length 
a) was observed, in good agreement with corollaries of the extended KH-a 
theorem and exact laws in MHD turbulence [22] • The LAMHD—a equa- 
tions also capture high-order statistics (up to and including order 8) and the 
anomalous scaling of the longitudinal structure function exponents, with a 
net gain in speed close to a factor of 16. For lower order structure functions, 
very little contamination of the scaling could be detected at scales larger than 
a. On the other hand, for the highest computed order, fluctuations in the 
scaling are observed for the runs with the smallest resolution. Note that we 
would not expect any scaling to be preserved for a so large that no inertial 
range remains. 

Turbulence closures are never unique. The present case may owe its suc- 
cess not to its particular form, but rather to its general properties that it 
(1) preserves physical avenues of nonlinear energy exchange and (2) allows 
correct vortex stretching. These two properties derive from its origin via a 
Lagrangian-averaged Hamilton's principle. The derivation also identifies the 
appropriate dissipation for proper energy decay, which involves an enhanced 
resistivity, but not an enhanced viscosity. Together, the Navier-Stokes viscos- 
ity and the enhanced resistivity produce regularization (e.g., existence and 
uniqueness of strong solutions and their global attractor of finite Hausdorf 
dimension, to be discussed elsewhere). In turn, these choices of viscosity and 
resistivity allow the intermittency found here, which might have otherwise 
been suppressed. 

Relying on the fact that, contrary to fluids, two dimensional MHD tur- 
bulence displays a direct cascade of energy and intermittency, we could show 
that the LAMHD—a equations reproduce intermittency features of turbu- 
lent flows and thus we postulate that these results will carry over to the 
three-dimensional case and thus these results could be also of relevance to 
the modeling of neutral fluids. Future challenges will include implementa- 
tion of the LAMHD—a model in domains with boundaries and the study 
of intermittency for magnetic Prandtl numbers besides unity. The choices 
of boundary conditions may be expected to strongly influence the solution 
behavior. Of course, this matter is beyond the scope of the present article. 
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